All-electron quantum Monte Carlo calculations for the noble gas atoms He to Xe 
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We report all-electron variational and diffusion quantum Monte Carlo (VMC and DMC) calcu- 
lations for the noble gas atoms He, Ne, Ar, Kr, and Xe. The calculations were performed using 
Slater- Jastrow wave functions with Hartree-Fock single-particle orbitals. The quality of both the op- 
timized Jastrow factors and the nodal surfaces of the wave functions declines with increasing atomic 
number, Z, but the DMC calculations are tractable and well-behaved in all cases. We discuss the 
scaling of the computational cost of DMC calculations with Z. 

PACS numbers: 02.70.Ss, 31.25.Eb, 71.15.Nc 



I. INTRODUCTION 

The variational quantum Monte Carlo (VMC) method 
and the more sophisticated diffusion quantum Monte 
Carlo (DMC) method^ can yield highly accurate ener- 
gies for many-electron systems. One of the main attrac- 
tions of these methods is that the cost of calculating the 
energy of N quantum particles scales roughly as N 2 - 
N 3 , which is better than other many-body wave function 
techniques. However, although the scaling with particle 
number is quite advantageous, the cost increases rapidly 
with the atomic number Z of the atoms involved. The- 
oretical estimates of this scaling 2,3 for DMC calculations 
have varied from Z 5 5 to Z 6 - 5 , while a practical tests in- 
dicated a scaling of about Z 5 - 2 . 

Numerous all-electron DMC studies have been re- 
ported^i&MO^^ 3 .^^^! for atoms up to Z = 10, 
but very few have included heavier atoms. DMC stud- 
ies of heavier atoms have normally used pseudopoten- 
tials to remove the chemically inert core electrons from 
the problem. However, pseudopotentials inevitably intro- 
duce some errors and it may be useful to consider how 
much progress can be made with all-electron DMC calcu- 
lations. Accurate all-electron calculations for atoms may 
also be useful in constructing pseudopotentials which in- 
corporate many-body effects. In this paper we report 
VMC and DMC calculations for the noble gas atoms He, 
Ne, Ar, Kr, and Xe, which extends the range of atoms 
studied within VMC and DMC up to Z = 54. The main 
aims of this paper are to investigate how well current all- 
electron DMC methods perform for heavy atoms and to 
study the scaling of the computational cost with Z . 



II. VMC AND DMC METHODS 

In VMC the energy is calculated as the expectation 
value of the Hamiltonian with an approximate many- 
body trial wave function containing a number of variable 
parameters. In DMC the estimate of the ground state 
energy is improved by performing an evolution of the 
wave function in imaginary timei The fermionic sym- 
metry is maintained by the fixed-node approximation 18 , 



in which the nodal surface of the wave function is con- 
strained to equal that of a trial wave function. Our DMC 
algorithm is essentially that of Umrigar et al£, and we 
employ the modifications to the Green function for all- 
electron calculations proposed in that paper. All of our 
VMC and DMC calculations were performed using the 
CASINO codej±2 

Our trial wave functions were of the standard Slater- 
Jastrow form, 



* = e DtD 



1^1 ■ 



(1) 



The Jastrow factors, e J , 
the variables r« = |r, 



were chosen to be functions of 
r, 3 — |±j — Tj| and ri = \vi\, where is 
the position of electron i with respect to the nucleus. 
Our Jastrow factors 20 for He, Ne, Ar, Kr, and Xe con- 
tained a total of 26, 75, 79, 80, and 54 adjustable pa- 
rameters, respectively^ 3 - The optimal parameter values 
were obtained by minimizing the variance of the energy 
within a VMC procedure! 21 : 22 The Slater determinants, 
D a , were formed from single-particle orbitals obtained 
from Hartree-Fock (HF) calculations using (i) numerical 
integration on a radial grid, and (ii) Gaussian basis sets 
and the CRYSTAL98 code. 24 Although the numerical or- 
bitals are the more accurate they are not available for 
molecular systems, in which Gaussian basis sets are very 
commonly used. 

In both VMC and DMC the energy is calculated as an 
average over a set of electron configurations of the local 
energy, El = ty~ 1 H'fy , where H is the Hamiltonian. The 
presence of core electrons causes two related problems. 
The first is that the shorter length scale variations in 
the wave function near a nucleus of large Z require the 
use of a small timestep. The second problem is that the 
fluctuations in the local energy tend to be large near the 
nucleus, because both the kinetic and potential energies 
are large. Although these fluctuations can be reduced by 
optimizing the trial wave function, in practice they are 
large for heavier atoms. 

At a nucleus the exact wave function has a cusp 2 ^ such 
that the divergence in the potential energy is canceled by 
an equal and opposite divergence in the kinetic energy. 
A determinant of exact HF orbitals obeys the electron- 
nucleus cusp condition. However, Gaussian functions are 
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smooth, and a determinant of such orbitals cannot have a 
cusp, so the local energy diverges at the nucleus. In prac- 
tice one finds wild oscillations in the local energy close to 
the nucleus, which increase the variance of the energy in 
VMC calculations and lead to timestep errors and even 
numerical instabilities in DMC calculations. To solve this 
problem we make small corrections to the single-particle 
orbitals close to the nucleus, which impose the correct 
cusp behavior. 26 



III. RESULTS 
A. Quality of the trial wave functions 

The main results of our HF, VMC, and DMC calcu- 
lations are shown in Table fl] Our HF energy for He 
with Gaussian orbitals is very close to the result with 
the "exact" numerical orbitals, indicating the high qual- 
ity of the Gaussian basis set used. For Ne the HF energy 
with Gaussian orbitals is a little higher than the value 
with the numerical orbitals, although this difference is 
not large enough to affect the DMC results. We experi- 
mented with various Gaussian basis sets for the heavier 
noble gas atoms and found that the basis set errors at 
the HF level tend to increase significantly with atomic 
number. In the case of Xe our best Gaussian basis set 
gave an error of 0.11 a.u. For Ar, Kr, and Xe we therefore 
used only the numerical orbitals. 

The exact ground-state wave function of a two-electron 
atom is a nodeless function of r\, r2, and r\2, which is the 
same form as our trial wave function for He. We therefore 
expect to obtain a highly accurate trial wave function 
for He. We refer to the difference in the HF and DMC 
energies as the "DMC correlation energy" . If one keeps 
the orbitals fixed and varies the Jastrow factor, then the 
lowest energy one could obtain is the DMC energy. The 
percentage of the DMC correlation energy retrieved at 
the VMC level is therefore a measure of the quality of 
the Jastrow factor. From the data in Table|T]we find that 
our VMC calculations retrieve 99.5%, 91%, 85%, 70%, 
and 59% of the DMC correlation energy for He, Ne, Ar, 
Kr, and Xe, respectively. We believe that the decrease 
in the quality of the Jastrow factor with increasing Z 
arises from the increasing inhomogencity of the atoms. 
Creating accurate Jastrow factors for all-electron studies 
of heavy atoms is a challenging problem. 

Our VMC and DMC energies for Ne obtained with 
the numerical orbitals are very close to those obtained 
in our earlier work-S. Huang et al^ obtained a VMC 
energy of —128.9008(1) a.u., which is only slightly lower 
than our value, although they also optimized the orbitals. 
Our DMC energies for Ne are within error bars of those 
reported by Umrigar et al£, but the remaining fixed-node 
error of 0.016 a.u. is substantial. 

From the data in Table Q] we observe that the percent- 
ages of the correlation energy missing at the DMC level 
are 0%, 4%, 9%, 18% and 23% for He, Ne, Ar, Kr and 



Atom 


Method 


Orb. 


Total energy 


E c 






type (a.u.) 






HF 


G 


-2.86165214 


% 




HF 


N 


-2.86168000 


% 




VMC 


G 


-2.903499(8) 


99.5 % 


He 


VMC 


N 


-2.903527(9) 


99.5 % 




DMC 


G 


-2.903732(5) 


100 % 




DMC 


N 


-2.903719(2) 


100 % 




"Exact"^ 




-2.903724 


100 % 




HF 


G 


-128.53832860 


% 




HF 


N 


-128.54709811 


% 




VMC 


G 


-128.8794(4) 


85 % 


Ne 


VMC 


N 


-128.891(5) 


88 % 




DMC 


G 


-128.9232(5) 


96 % 




DMC 


N 


-128.9231(1) 


96 % 




"Exact"^ 


. 


-128.939 


100 % 




HF 


N 


-526.81751277 


% 


Ar 


VMC 


N 


-527.3817(2) 


77 % 




DMC 


N 


-527.4840(2) 


91 % 




"Exact"-22 




-527.55 


100 % 




HF 


N 


-2752.05497715 


% 


Kr 


VMC 


N 


-2753.2436(6) 


57 % 




DMC 


N 


-2753.7427(6) 


82 % 




"Exact"^ 




-2754.13 


100 % 




HF 


N 


-7232.13836331 


% 


Xe 


VMC 


N 


-7233.700(2) 


46 % 




DMC 


N 


-7234.785(1) 


77 % 




"Exact"-^ 




-7235.57 


100 % 



TABLE I: Total energies of the noble gas atoms, and the 
percentages of the correlation energies, E c , retrieved. (G) de- 
notes a calculation with a Gaussian basis set and (N) denotes 
numerical orbitals. The "exact" energies were obtained from 
data in the indicated references. 

Xe, respectively. This indicates that the size of the fixed- 
node error increases rapidly with Z. 

B. Theoretical scaling with atomic number 

It is of interest to study the CPU time required to 
obtain a fixed standard error in the mean energy, A, as 
a function of the atomic number Z. The required CPU 
time, T, can be written as 

TocMTc, (2) 

where M is the total number of generations of electron 
configurations and Tq is the CPU time for one move of 
C configurations, where C is the average number of con- 
figurations in a generation. 

Ceperley 3 showed that A 2 can be written as the sum 
of two terms; the first corresponds to the square of the 
standard error evaluated as if the DMC energies were 



uncorrelated, and the second accounts for the effects of 
correlations. In DMC calculations the timestep r is nor- 
mally chosen to be small, and the correlations between 
configurations at successive generations are large, so that 
the second of these terms dominates. Ceperley showed 
that this term is given approximately by 3 - 



2|£Vmc — £dmc| 
tMC 



(3) 



Since A 2 is inversely proportional to the total number of 
configurations, we obtain 



T oc t~ l \E- 



VMC 



-E 



DMC 



T c ■ 



(4) 



Ceperley used the simple approximation -Evmc — 
-EdmcI E c , and the approximate scaling E c oc Z 15 . 
He also argued that avoiding large timestep errors re- 
quires r oc Z~ 2 , as the average distance diffused should 
be smaller than the size of the Is orbital, which is pro- 
portional to Z^ 1 . Finally, he used Tc oc Z 2 to obtain an 
overall scaling of 



T oc Z 



5.5 



(5) 



Hammond et al. 2 argued along similar lines, although 
they chose Tc oc Z 3 , leading to an overall scaling of T oc 
Z 65 , In what follows we examine some aspects of these 
arguments. 



C. Numerical tests of scaling with atomic number 

Discussions of the actual scaling of the computational 
cost of calculations with system size or atomic number 
are fraught with difficulties. The results depend on the 
computers on which the calculations are run, the algo- 
rithms used, and on the details of the software used. Our 
calculations are run on parallel computers in which each 
processor deals with a small number of electronic config- 
urations (one in VMC and roughly ten in DMC). The in- 
terprocessor communications are negligible in VMC and 
small in DMC, and the computational cost is inversely 
proportional to the number of processors used. All of 
the DMC results used for determining the scaling of the 
computational cost with atomic number were performed 
on 96 processors of a Sunfire Galaxy machine, although 
most of the variance minimizations were performed on a 
cluster of 16 Xeon dual-processors. 

To ensure that timestep errors are small the DMC 
timestep should be chosen so that the probability of a 
move being accepted is high. For the DMC results re- 
ported in Table U we used timesteps of 0.02, 0.0025, 
0.0009, 0.00035, and 0.0002 a.u., for He, Ne, Ar, Kr, 
and Xe, respectively, which were chosen so that in each 
case slightly more than 99 % of the proposed moves were 
accepted. These timesteps scale as Z~ 1A1 , which is sig- 
nificantly weaker than the Z~ 2 scaling used in the earlier 
theoretical estimates We therefore expect the timestep 
bias in our DMC results to increase with Z. 
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FIG. 1: The correlation energy E c as a function of atomic 
number Z. Crosses: "exact" values, diamonds: DMC values. 
The dotted line is a fit to the "exact" values giving E c oc 
Z s , while the dashed line is a fit to the DMC data giving 
E c oc Z 1 - 26 . 



We thoroughly investigated the timestep dependence 
of the energies for He and Ne, concluding that they are 
negligible compared with the statistical error bars given 
in Tabled For each of Ar, Kr and Xe we performed calcu- 
lations at four different timesteps, and we estimate that 
the timestep errors in the corresponding DMC energies 
are less than 0.002 a.u. (Ar), 0.01 a.u. (Kr) and 0.015 a.u. 
(Xe). It is likely that the larger timestep errors in our 
DMC results for the heavier atoms arise both from the 
reduction in the quality of the trial wave functions and 
the poorer sampling of the core electrons. 

The correlation energy is normally defined as the dif- 
ference between the exact non-relativistic ground state 
energy and the HF energy, assuming static point nuclei. 
Accurate estimates of the correlation energies of neutral 
atoms for Z — 2 to 18 are given by Davidson and Chakra- 
vorty 2 ^, while Clementi and Hofmanr* 3 -^ give values for Kr 
and Xe which, while probably not as accurate as those 
for the lighter atoms, are expected to be quite reliable. 
We will take these as our reference data and refer to them 
as the "exact" correlation energies, and when added to 
the Hartree-Fock energies, the "exact" energies. 

Fig- [H shows the correlation energy as a function of Z 
from our DMC data and the estimates of Davidson and 
Chakravorty 2 ^ (Z=2-18), and Clementi and Hofmann 30 
(Z = 36,54). It is clear that (apart from He) DMC 
underestimates the correlation energy, and that the un- 
derestimation becomes more severe at larger Z. The best 
power-law fit to the "exact" data for the noble gas atoms 
gives E c oc Z 1,33 , while for our DMC data we obtain 
Z 126 . The scaling of Z 1 - 5 assumed in the earlier theo- 
retical estimates^ 3 - is somewhat of an overestimate. 

As mentioned in Sec. IIII B[ the quantity which actu- 
ally enters Ceperley's approximation of Eq. ^ for the 
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Atom |rMCA 2 (a.u.) |_E VM c - £dmc| (a.u.) 



He 0.00019 0.00019 

Ne 0.024 0.032 

Ar 0.11 0.10 

Kr 0.66 0.50 

Xe 1.3 1.0 



TABLE II: The quantity \tMCA 2 and the difference be- 
tween the VMC and DMC energies. 

20i ■ 1 ■ 1 ■ 1 ■ 1 1 




■ | , 1 , 1 , L 

12 3 

ln(Z) 



FIG. 2: The logarithm of the CPU time required to obtain 
a fixed error bar in the energy versus ln(Z) for our DMC 
calculations. The dashed line shows the fitted scaling of Z 5 47 . 
The CPU times are measured in seconds. 

variance of the DMC energy is the difference between 
the variational and DMC energies, |-EVmc — ^dmcI- Us- 
ing the VMC and DMC data given in Table [fl we find 
I-EVmc — -EdmcI oc Z 262 . The reason that this quantity 
increases more rapidly with Z than E c is that the per- 
centage of the correlation energy retrieved in our VMC 
calculations decreases with Z more rapidly than in our 
DMC calculations. 

We can also test Eq. Q directly by comparing the 
difference between the variational and DMC energies 
I-Bvmc — -EdmcI with the quantity \tMCA 2 . Serial cor- 
relation of the data has been taken into account when 
computing the variance over the run. The results shown 
in Table [IT] indicate that the two quantities are in good 
agreement, which is rather satisfactory considering the 
large range of Z and the very different qualities of trial 
wave functions used. The quantity ^tMCA 2 is fitted by 
a scaling of Z 2 - 71 . 

We found that the computational cost of moving all the 
electrons in a configuration scaled as Z 135 in our DMC 
calculations. This is rather better than the scalings as- 
sumed by Ceperley^ {Z 2 ) and by Hammond el al& {Z 3 ). 
If we studied a system containing many atoms, the scal- 
ing of the computational cost for moving all the electrons 
in a configuration would be expected to increase roughly 



as N 2 , although the use of localized Wannier functions 
could reduce this to NM 

Putting together our scalings for the factors in Eq. (|4]) 
we find 

T oc Z 1A1 x Z 2m x Z 135 = Z 5 - 3S . (6) 

We can now compare this with the actual DMC com- 
putations reported in Table Q] In Fig. [2] we show the 
logarithm of the CPU time as a function of Z for a given 
standard error of the mean. The best fit gives a scaling 
of Z 5A7 , in good agreement with the prediction of Z 538 
from Eq. ([4} . 

As mentioned before, our DMC results for the heavier 
atoms suffer from significant timestep errors. If we adopt 
the Z~ 2 scaling for the timestep instead of the Z~ 1A1 
used above, we obtain an overall scaling of T oc Z 5 97 , 
which is higher than the value of T oc Z 5 2 obtained in 
the practical tests of Hammond et alA Moreover, it seems 
likely that an even more rapid scaling would be required 
to achieve a timestep error independent of Z . 



IV. CONCLUSIONS 

We have applied the VMC and DMC methods to no- 
ble gas atoms up to Xe (Z = 54), using Slater- Jastrow 
wave functions with Hartree-Fock single-particle orbitals. 
The percentage of the DMC correlation energy obtained 
at the VMC level decreases with Z, indicating that the 
quality of our Jastrow factors decreases with Z. The 
percentage of the exact correlation energy retrieved at 
the DMC level also decreases with Z, indicating that the 
quality of the HF nodal surface deteriorates with increas- 
ing Z. 

Our study shows that Ceperley's expression^ for the 
variance of the DMC energy (Eq. ([3])) is accurate to bet- 
ter than a factor of two for the systems studied here. The 
computational cost required to obtain a fixed statistical 
error bar in the energy scaled as Z 5A7 , but in these calcu- 
lations the timestep error increased significantly with Z. 
The scaling required to achieve a timestep error indepen- 
dent of Z is difficult to estimate, but it would certainly 
be higher than Z 5AJ . However, it may well be reasonable 
to incur substantial timestep errors deep in the core of 
the atom when we calculate chemical properties which 
are related to the valence electrons. 
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